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ABSTRACT 

The observed sample of double neutron-star (NS-NS) binaries presents a challenge to population- 
synthesis models of compact object formation: the input model parameters must be carefully chosen 
so the results match (i) the observed star formation rate and (ii) the formation rate of NS-NS binaries, 
which can be estimated from the observed sample and the selection effects related to the discover- 
ies with radio-pulsar surveys. In this paper, we select from an extremely broad family of possible 
population synthesis models those few (2%) which are consistent with the rate implications of the ob- 
served sample of NS-NS binaries. To further sharpen the constraints the observed NS-NS population 
places upon our understanding of compact-object formation processes, we separate the observed NS- 
NS population into two channels: (i) merging NS-NS binaries, which will inspiral and merge through 
the action of gravitational waves within 10 Gyr, and (ii) wide NS-NS binaries, consisting of all the 
rest. With the subset of astrophysically consistent models, we explore the implications for the rates 
at which double black hole (BH-BH), black hole-neutron star (BH-NS), and NS-NS binaries merge 
through the emission of gravitational waves. 

Subject headings: binaries: close — stars:evolution — stars:neutron - black hole physics - stars:winds 



1. INTRODUCTION 

Interest in the formation channels and rates of dou- 
ble compact objects (DCOs) has increased in recent 
years partly because, at the late stages of their inspi- 
ral through the emission of gravitational waves, they can 
be strong enough sources to be detected by the many 
presently-operating ground-based gravitational-wave de- 
tectors (i.e., LIGO, GEO, TAMA). But with the notable 
exception of NS-NS mergers - see Kim, Kalogera, and 
Lorimer (2003), henceforth denoted KKL - the merger 
rates for DCOs with black holes have not been con- 
strained empirically. The only route to rate estimates 
for double black hole (BH-BH) and black hole-neutron 
star (BH-NS) binaries is through population synthesis 
models. These involve a Monte Carlo exploration of the 
likely life histories of binary stars, given statistics gov- 
erning the initial conditions for binaries and a method 
for followi ng the behavior o f single and binary stars (see, 
e.g.. iBelczvnski et al.ll20 021 Unfortunately, our under- 
standing of the evolution of single and binary stars is 
incomplete, and we parameterize that uncertainty with 
a great many parameters (~ 30), many of which can 
cause the predicted DCO merger rates to vary by more 
than an order of magnitude when varied independently 
through their plausible range. To arrive at more defini- 
tive answers for DCO merger rates, we must substan- 
tially reduce our uncertainty in the parameters that enter 
into population synthesis calculations through compari- 
son with observations. 

Clues to the physics underlying the formation of tight 
compact binaries can be obtained through a study of 
each individual known DCO system. Some authors have 
followed this path, for example examining the potential 
evolutionary and kinematic histories of each individual 
binary to deduce the pulsar kicks needed to rep roduce 
their evolutionary path (e.g., Willcm s et a.lJl2f)Q4jfl . How- 
ever, the simplest and most direct way to constrain the 
parameters of a given population synthesis code is to 
compare several of its many predictions against observa- 



tions. For example, the empirically estimated formation 
rates derived from the six known Galactic NS-NS bina- 
ries - half of which are tight enough to merge through the 
emission of gravitational waves within 10 Gyr - should be 
reproduced by any physically reasonable combination of 
model parameters for population synthesis. 

In this paper, we describe the constraints the observed 
NS-NS population places upon the most significant pa- 
rameters t hat enter into one population synthesis code, 
StarTrack l|Belczvnski et alJ 120021 120051) . Furthermore, 
we use the set of models consistent with this observa- 
tional constraint to revise our population-synthesis-based 
expectations for various DCO merger rates. 

In Sectional we describe the observational constraints 
from NS-NS: reviewing and extending the work of Kim, 
Kalogera, and Lorimer (2003) we briefly summarize the 
observed sample of NS-NS binaries, the surveys which 
detected them, and the implications of the (known) sur- 
vey selection effects for the expected NS-NS formation 
rate. In Section we describe our population synthe- 
sis models and their predictions for NS-NS binary (and 
other compact binary) formation rates. Since a compre- 
hensive population synthesis survey of all possible mod- 
els is not computationally feasible, we describe an effi- 
cient approximate fitting technique [used previously in 
O'Shaughnessy, Kalogera, and Belzcynski (2005); here- 
after OKB] we developed to accurately approximate the 
results of complete population synthesis calculations. Fi- 
nally, in Section 21 we select from our family of possible 
models those predictions which are consistent with the 
observational constraints. We then employ that sample 
to generate refined predictions for the expected BH-BH, 
BH-NS, and NS-NS merger rate through the emission of 
gravitational waves. 

We find the observed NS-NS population can provide a 
tight constraint (albeit a complicated one to interpret) 
on the many parameters entering into population synthe- 
sis models. With this article serving as an outline of the 
general method, we propose to impose in the future sev- 
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eral additional constraints, including notably the lack of 
any observed BH-NS systems, the empirical supernova 
rates as well as formation rates of binary pulsars with 
white dwarf companions. 

2. EMPIRICAL RATE CONSTRAINTS FROM THE 
NS-NS GALACTIC SAMPLE 

Seven NS-NS binaries have been discovered so far in 
the Galactic disk. Recently KKL developed a statisti- 
cal method to calculate the probability distribution of 
rate estimates derived using the observed sample and 
modeling of survey selection effects. Four of the known 
systems will have merged within lOGyr (i.e., "merging" 
binaries: PSRs J0737-3039, B1913+16, B1534+12, and 
J1756-2251) and three are wide with much longer merger 
times (PSRs J1811-1736, J1518+4904, an d J1829+2456). 
PSR J1756-2251 was discovered recently ((Faulkner et all 
12005(1 and has not been included in the current calcula- 
tions. We do not expect, however, that this system will 
sign ificantly change our exp ectations of the merger rate. 
(see lKalogera et all l|2004aD for details). This fourth sys- 
tem is sufficiently similar to PSR B1913+16 and was dis- 
covered with pulsar acceleration searches, the selection 
effects of which have already been accounted for (see Ta- 
ble for the properties of the six systems used here) . 

In what follows we use observational constraints 
based on the r ate prob ability distribution derived by 
iKalogera et alJ l)2004alb[) for merging binaries and the 
equivalent results for the three wide binaries (presented 
for the first time in this article). In what follows we use 
the index k = 1 . . . 6 to refer to the three merging (1,2,3) 
and the three wide (4,5,6) NS-NS binaries. 

2.1. Merging NS-NS Binaries 

The discovery of NS-NS binaries with radio pulsar 
searches and our understanding of the selection effects 
involved allows us to estimate the total number of such 
systems in our Galaxy and their formation rate. KKL 
developed a statistical analysis designed to account for 
the small number of known systems and associated un- 
certainties. Specifically, they found that the posterior 
probability distribution function Vk for NS-NS formation 
rates IZk for each sub-population k of pulsars similar to 
the fcth known binary pulsar is given by 

v k {n) = Alne- A * n . (i) 

The parameter Ak depends on some of the properties 
of the pulsars in the observed NS-NS sample [see KKL 
Eq. (17)]: 

A = Tiife/ (fbNpsn) (2) 

where is the fraction of all solid angle the pulsar 
beam subtends; Tuf e is the total binary pulsar lifetime 

Tufe = r sd + T rnrg (merging) (3) 

(where r st j is the pulsar spindown age [see Arzoumanian, 
Cordes, & Wasserman 1999] and T mrg is the time re- 
maining until the pulsar merges through the emission 
of gravitational waves [see Peters 1964 and Peters and 
Mathews 1963]); and Npsr is the total estimated num- 
ber of systems similar to each of the observed one (i.e., 
NpsR * s effectively a volume-weighted probability that 
a pulsar with the same orbit and an optimally oriented 



beam would be seen with a conventional survey; this fac- 
tor incorporates all our knowledge of pulsar survey se- 
lection effects as well as the pulsar space and luminosity 
distributions). Tabled lists f° r each merging NS-NS bi- 
nary several intrinsic parameters (i.e., the best known 
values for /;,; several lifetime-related parameters, such 
as r s d and T mgr ) and two key quantities which depend 
on our analysis of selection effects: Npsr and the de- 
duced A [i.e., via Eq. J3J]. [Results are shown for our 
preferred model for binary pulsar space and luminosity 
distribution; sec model #6 and details in KKL.] The 
total NS-NS posterior density of the combined rate rep- 
resented by the observed samples can be computed by a 
straightforward convolution, 

v{n tot ) = J dn 1 dn 3 d'R25{ntot -Hi-Hi- n 3 ) 

x Vx (^1)^2(^2)^3 (fts) (4) 

described in detail in Section 5.2 of KKL and p resented in 
detail f or the three-binary case in Eq. (A8) of lKim et aU 
l!2004bH . 

KKL also demonstrated that the resulting rate distri- 
butions depend only weakly on the spatial distribution of 
NS-NS locations (see their Figure 7 and the end of their 
Section 6). Thus the NS-NS rate distribution effectively 
depends on only one model assumption, the choice of the 
intrinsic radio pulsar luminosity function - which, in the 
KKL approach is given by [see KKL Eq. (3), following 
Cordes and Chernoff (1997)], 

4>{L)dL = (p- l)(L/L min )-PdL/L min . (5) 

Thus it is controlled by two parameters, the minimum 
allowed pulsar luminosity [L m in) and the power law p > 
1 governing their relative luminosity probabilities. 

KKL did not complete their calculation for a compre- 
hensive posterior probability distribution for the NS-NS 
rate estimates, however, because up-to-date empirical 
pro bability constrai nts for p and L m i n are not available 
Ccf . . IKalogera! (|2004T) ) . Instead, they presented results for 
a few selected models, emphasizing one model (model 6) 
whose properties (L m in = 0.3mJy (kpc) 2 and p = 2) are 
close to the median values they expect will be found when 
all present observations are taken into account. For this 
particular model, the empirical parameters Ak which de- 
scribe the posterior densities are given in Tabled 

2.2. Wide NS-NS Binaries 

The same general technique outlined above can be ap- 
plied to the formation rate of wide NS-NS binaries: the 
same form of distribution function V(TZ) [Eq. G] applies 
and it depends on the same parameter A [Eq. |2] . The 
main change is the relevant lifetime. Since these binaries 
do not merge, their detectable lifetime is now the sum 
of the time remaining before the pulsar spins down (r s d, 
described earlier) and the length of time the pulsar will 
remain visible (the "death time" of a pulsar; see Chen & 
Ruderman 1993). However, since T s d estimates are some- 
what uncertain, we require that they do not exceed the 
current age of the Galactic disk (lOGyr). To summarize, 
then, the only change from the previoius approach is to 
replace the previous expression for the lifetime, Eq. ©, 
with 

Tufe = min (r a d, lOGyr) + r d (wide) . (6) 
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PSRs 


pa 

s 




(ms) 


(1) merging NS-NS 


B1913+16 


59.03 


B1534+12 


37.90 


J0737-3039 


22.70 


(2) wide NS-NS 


J1811-1736 


104.182 


J1518+4904 


40.935 


J1829+2456 


41.0098 



a Spin period. 
b Spin-down rate. 



c Orbital period 

d Estimated mass of a companion, where Mns is assumed to be 1.35 Mq 
e Eccentricity. 

f Characteristic age of a pulsar. 

g Spin-down age of a pulsar. We calculate t s( j only for those recycled pulsar. \ 5 
^Merging time of a binary system due to the emission of gravitational waves. 
'Death time of a pulsar. | 
J Most probable value of the total number of pulsars in a model galaxy estimated fop the refer e 
k Beaming factor for pulsar Q ^ 

'Parameter in rate equation [see Eq. Gil. 
""References: (1) Lundgren, Zepka, & Cordes (1995); (2) Edwards, & Bailes (200J) 
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|t (2003) ; (5) van Kerkwij] 



Table ^ lists pulsar parameters and deduced quantities 
for the three wide NS-NS binaries used in this study. 
Current pulsar observations do not provide us with any 
estimates of the beaming fractions relevant to the pulsars 
in these wide systems. Guided by the beaming fraction 
distribution for merging pulsars, we adopt a value of 6 for 
the beaming factor for all the wide NS-NS pulsars. For 
simplicity, we present the results for only the preferred 
luminosity model (i.e., for the specific choice for p and 
Lmin mentioned above). These distributions again follow 
Eq. Q , with parameters given by Table ^ where 
Ak is determined for each pulsar class k from physical 
parameters presented in the table. 

For each class s eparately (merging and wide binaries) 
we use Eq. (A8) of lKim et alJ l)2004b|) to generate a com- 
posite probability distributions for the formation rate 
of binaries in that class: V m (merging) and V w (wide). 
Thus we arrive at the two estimates shown in Fig. ^ f° r 
the empirical probability distribution 

p(logft) =V{n)Tl\nlQ 

for the formation rates 7Z m and 1Z W of these two classes 
of binary. ^From these distributions we derive a 95% 
confidence intervals for each formation rate; for example, 
the upper and lower rate limits 1Z W ,± satisfy 

I dUV w (n)= / dTZV w (TZ) = 0.025 . (7) 
Jo Jt^w,+ 
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Fig. 1. — A plot of empirically-deduced probability distributions 
for merging (right) and wide (left) NS-NS binaries; see Sec. [2] The 
solid vertical lines are at (i) log 10 U = -4.5388, -3.49477, the 95% 
confidence interval for the merging NS-NS merger rate; and at (ii) 
\og 1Q K = -6.7992,-5.7313, the 95% confidence interval for the 
wide NS-NS formation rate. 



3. ESTIMATES FOR MERGER RATES 
3.1. Population Synthesis Estimates 

We estimate formation and merger rates for several 
classes of double compact objects using the Star Track 
code first developed by Belczynski, Kalogera, and Bu- 
lik (2002) [hereafter BKB] and recently significantly up- 
dated and tested as described in detail in Belczynski et al. 
2005. In this code, seven parameters strongly influence 
compact object merger rates: the supernova kick distri- 
bution (3 parameters), the massive stellar wind strength 
(1), the common-envelope energy transfer efficiency (1), 
the fraction of mass accreted by the accretor in phases of 
non-conservative mass transfer (1), and the binary mass 
ratio distribution described by a negative power-law in- 
dex (1). To allow for an extremely broad range of possi- 
ble models, we used the specific parameter r anges quoted 
in Section 2 of IQ'Shaughnessv et al.1 l|2005() . 

We randomly choose model parameters in this space 
and evaluate their implications, by progressively exam- 
ining the evolution of binary after binary. We then ex- 
tract from our simulations predictions for several DCO 
formation rates (BH-BH, BH-NS, and NS-NS) by scaling 
up the ratio of DCO formation events we obtain in each 
simulation (n) to the total number of binaries studied in 
the simulation (N) by a factor proportional to the ex- 
pected ratio between N and the number of stars formed 
in the Milky Way. We set this scaling factor by assum- 
ing a constant star- formation rate of M m S.^Mqyt^ 1 , 
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as described in the Appendix of OKB. 

Extracting predictions for the "visible" NS-NS forma- 
tion rates: To compare the predictions of population syn- 
thesis calculations against the empirical rate constraints 
derived for the pulsar samples, we must determine the 
formation rates of NS-NS binaries that could be "vis- 
ible" as pulsars. Since we do not follow the detailed 
pulsar evolution with StarTrack (due to major uncer- 
tainties related to pulsar magnetic field evolution), we 
choose a minimal criterion for identifying NS-NS bina- 
ries that possibly contain a recycled pulsar: if the first 
NS in the binary has experienced any accretion episode 
(through either Roche-lobe overflow and disk accretion or 
a common-envelope phase), then the binary is identified 
as a potential binary recycled pulsar and is included in 
the calculation of the NS-NS "visible" pulsar formation 
rate. 

Practical complications in merger rate calculations: 
We would have preferred to proceed as in OKB and per- 
form, for each separate DCO type (e.g., BH-BH binaries), 
a sequence of Monte Carlo computations tailored to de- 
termine this type's merger rate to some fixed accuracy 
(say, 30%) as a function of all population synthesis pa- 
rameters. Instead, owing to computational limitations, 
we had to extract multiple types of information from 
each population synthesis run; Appendix lAl describes in 
greater detail the collection of population synthesis runs 
we performed and the manner in which these runs were 
used to estimate various DCO formation rates. 

3.2. Mapping population synthesis rates versus 
parameters 

In order to constrain population synthesis parameters 
based on rate measurements, we must be able to invert 
the relation between rate and model parameters to find 
all possible models consistent with a given rate. In other 
words, we must fit the rates over all seven parameters. 

OKB first demonstrated that, even using sparse data in 
a high-dimensional space of population synthesis param- 
eters, an effective fit could be found for formation rates 
of DCOs (see OKB Fig. 2 and their Section 4). We con- 
structed separate polynomial least-squares fits to each of 
the five rate functions we need (i.e., for BH-BH, BH-NS, 
visible merging NS-NS, and visible wide NS-NS binaries 
we performed a cubic least-squares fit; and for the overall 
NS-NS merger rate - including all merging NS-NS bina- 
ries, whether we expect them to be electromagnetically 
visible or not - we used a quartic least squares fit). Fig- 
ure |21 demonstrates that the fit is good: the errors are on 
the limiting scale we would expect, given the uncertain- 
ties in the input [i.e., the standard deviation of the loga- 
rithmic rate errors, \og w lZfi t /Htrue are 0.167 (BH-BH), 
0.22 (BH-NS), and 0.086 (NS-NS), are comparable with 
the minimum possible uncertainty we would expect given 
a perfect fit, log 10 (l + 1/Vl0) « 0.119; see Appendix lAl 
for a detailed discussion of the minimal uncertainties ex- 
pected for each rate] . The fits are sufficiently good that 
for our purposes we can replace population synthesis cal- 
culations with evaluations of our fits. In particular and 
by way of example, in Figure we generate a histogram 
for various DCO formation rates predicted by popula- 
tion synthesis, using (i) the actual outputs deduced from 
the population synthesis code, sampled at random Monte 
Carlo points (dashed line) , and (ii) the outputs obtained 



-8 -7 -6 -5 -4 




log 10 (R fit yr) 

Fig. 2. — This plot of the log 10 of the Galactic rate versus 
our fit to the rate, shown for BH-BH, BH-NS, and NS-NS sample 
points (all superimposed). The shaded region is offset by a factor 
1 ± 1/ ^/W. This region estimates the error expected due to random 
fluctuations in the number of binary merger events seen in a given 
sample. (See the appendix for a discussion of the number of sample 
points actually present in various runs.) 
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Fig. 3. — A plot of the a priori probability distribution for the 
BH-BH (left), BH-NS (center), and NS-NS (right) merger rates, 
versus the log 10 of the rate. These distributions were generated 
from the population synthesis code (dashed line) and fits (solid 
lines) assuming all parameters in the population synthesis code 
were chosen at random in the allowed region. 



from fits to the dataset from (i), sampled at a much larger 
number of data points (to insure smoothness; solid line). 
The two methods produce strikingly similar histograms, 
demonstrating that the fit will be adequate for our pur- 
poses. 

4. CONSTRAINTS FROM NS-NS OBSERVATIONS 

In Sec.|2]we constructed two straightforward empirical 
constraints (i.e., confidence intervals for the formation 
rate of "visible" merging and wide NS-NS binaries) we 
could place on the output of population synthesis calcu- 
lations. In this section, we apply these two constraints, 
individually and together, and determine their effect on 
rate predictions of population synthesis calculations. 

Bounding merging NS-NS rate: Figure^superimposes 
the observational bounds taken from the observed merg- 
ing NS-NS distribution (i.e., the 95% confidence interval; 
see Fig. ^| on top of the distribution of "visible" NS- 
NS merger rates obtained from unconstrained population 
synthesis. We limit attention only to those population 
synthesis models consistent with our constraint; specifi- 
cally, we randomly choose population synthesis models, 
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Fig. 4. — A plot of the a priori probability distributions for the 
visible merging (top) and visible wide (bottom) NS-NS formation 
rates produced from population synthesis. The dashed curves de- 
note the direct result from population synthesis; the solid curves 
denotes the result deduced from artificial data generated from a 
multidimensional fit to the visible wide and merging NS-NS rate 
data. The vertical lines are the respective 95% CI bounds presented 
in Fig. [3 



evaluate the "visible" merging NS-NS rate using our fit, 
and retain the model only if the rate lies within these two 
bounds. We find we reject 72% of models we initially 
considered plausible. For each of the small residual of 
consistent models, we can evaluate the BH-BH, BH-NS, 
and NS-NS merger rates (again using our fits). We find 
the merger rates increase slightly on average: the mean 
merger rate increases by a factor xl.2 for BH-BH, x2 
for BH-NS, and x2.2 for NS-NS. 

Bounding wide NS-NS rate: Figure 21 also shows the 
95% confidence interval for the wide visible NS-NS for- 
mation rate (Fig. ^) on top of our a priori population 
synthesis distribution for the wide visible NS-NS forma- 
tion rate. We find that with this constraint we have 
to exclude 80% of a priori plausible models. Using only 
models which satisfy this second constraint, we find DCO 
merger rates have dropped relative to our a priori predic- 
tions: the average BH-BH, BH-NS, and NS-NS merger 
rates are reduced by a factor 0.65, 0.26, and 0.22, respec- 
tively. 

Both constraints simultaneously: Very few population 
synthesis models (less than 2%) satisfy both constraints 
simultaneously. Since the set of consistent models is 
much smaller than initially permitted a priori, we have 
less uncertainty in our predictions: the standard devia- 
tion in the log of the merger rates has changed from our 
initial a priori uncertainty of 0.60, 0.55, and 0.63 (i.e., 
plus or minus a factor of 4, 3.5, and 4.3) for BH-BH, NS- 
NS, and BH-NS mergers, respectively, to 0.61, 0.44, and 
0.48 (i.e., plus or minus a factor of 4, 2.8, and 3). Further, 
since the wide constraint proves slightly more restrictive, 
the mean merger rates have dropped slightly from our 
prior expectations: the average predicted merger rates 
are 1.1/Myr [BH-BH] (down by a factor 0.43), 1.4/Myr 
[BH-NS] (down by a factor 0.24), and 6.7/Myr [NS-NS] 
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Fig. 5. — A plot of the a priori probability distribution for 
the BH-BH (top), BH-NS (center), and NS-NS (bottom) merger 
rates per milky way equivalent galaxy. As in Fig. [3] the dashed 
curves show the results obtained from our population synthesis 
calculations (i.e., our raw code results, smoothed); the thick solid 
curves show the results after we impose both our observational 
constraints (i.e., consistency with the observed number of visible 
wide and visible merging NS-NS binaries). 

(down by a factor 0.25). 

4.1. Advanced LIGO detection rates 

While we have presented the number of mergers occur- 
ring per Milky Way equivalent galaxy, advanced LIGO's 
inspiral detection range depends on the masses of the 
component objects. Specifically, one 4 km advanced 
LIGO detector (see Harry 2005) is expected to detect a 
binary with chirp mass M c = (TO 1 m 2 ) 3 ^ 5 /(rni + m^) 1 / 5 
(at signal-to-noise ratio 8) out to a distance 

d = 191 Mpc (M C /1.2M Q ) 5/6 . (8) 

Thus, if the real merger rate and chirp mass distribution 
for type a are lZ a and p a (M c ), respectively, then LIGO 
will on average detect a merger events at a rate 

ftaxico = 0.042 K a ((.M C /M Q ) 15 / 6 ) Q (9) 

where ((M C ) 15/6 ) Q = / dM cPa (M c )Ml 5/6 , and where 
for simplicity we assume a uniform distribution of Milky 
Way equivalent galax ies with density 0.01/(Mpc) 3 (see 
iNutzman et alJ l)2004l) for a discussion of short-scale cor- 
rections to this distribution in the case of short-range 
interferometers, like initial LIGO). 

Since Fig. El was produced using merger rate fits, we 
do not have the chirp mass information needed to trans- 
late that figure into a corresponding distribution for 
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Fig. 6. — This plot shows our expectations for LIGO's detection 
rates for merging BH-BH (solid line), BH-NS (dashed line), and 
NS-NS (dotted line) binaries. This plot was obtained directly from 
Fig. fusing Eq. J^j. 



the LIGO detection rate. In what follows we adopt 
mean chirp masses as derived from all our models in the 
archives. We find: 

( ( ^ c)15/6 )bh-bh = 111M ® /6 (10) 
((.M c ) 15/6 \ =5.8a4 5/6 (11) 

\ /BH-NS U 

((M c ) 15 / 6 ) =2Mg /6 . (12) 

\ /NS-NS 

These are to be compared to 224Mq 5 ^ 6 for two 10 M© 

BH, to 15.5Af 5/6 for a 10 M Q BH and a 1.4 M Q NS, and 
to 1.2M Q for two 1.4M Q NS. Figure presents our pre- 
liminary estimates for the advanced LIGO detection rate 
distribution. Note this figure provides the same informa- 
tion as Figure 03 except that the merger rates for each 
species have been rescaled according to Eq. ©. 

5. SUMMARY AND CONCLUSIONS 

In the context of matching theory with observations of 
the binary NS population, we have described how to con- 
strain the predictions for DCO merger rates population 
synthesis codes such as StarTrack by using two specific 
observational constraints. We find that to be consistent 
with the rate statistics of the observed NS-NS population 
(at the 95% confidence interval) , we must exclude at least 
98% of all the models we think a priori likely. We do not 
focus on explicitly describing the seven-dimensional re- 
gion of StarTrack model parameters consistent with our 
constraints, both because (i) we lack a compact way to 
describe a seven-dimensional region, and (ii) our region 
has meaning specifically for the StarTrack code. Other 
codes have different parameterizations of the same physi- 
cal phenomena, leading to potentially quite different rep- 
resentations of the same constraint region. However, in 
Appendix [0 we give some information about the mean 
constraints on population synthesis parameters but also 
the strong variance around these mean values. As de- 
scribed in Section 01 and particularly via Fig. to ex- 
tract a physically meaningful statement about the effect 
of imposing constraints (as opposed to a describing in- 
formation about parameters of one particular code), we 



have described how these constraints have improved our 
understanding of three DCO merger rates [BH-BH, BH- 
NS, and NS-NS]. We find that, using these two initial 
constraints, (i) the most probable merger rates (i.e., at 
peak probability density) are systematically lower than 
we would expect a priori, at least by a factor 2; and (ii) we 
reduce the uncertainty in the BH-NS and NS-NS merger 
rates by moderate factors (i.e., the standard deviation of 
log 1Z drops by 0.03 and 0.16, respectively). 

This paper only outlines the beginning of a large pro- 
gram we have undertaken to better constrain our under- 
standing of the evolution of single and binary stars and 
the associated predictions for gravitational-wave sources. 
We intend to add a few additional empirical constraints 
of DCOs (e.g., WD-NS binaries) and the lack of obser- 
vations of certain binary compact objects (notably, BH- 
NS binaries). Apart from rate constraints, the observed 
properties of DCOs (mass ratios, orbital separations and 
eccentricities) could also be used as constraints. Fur- 
ther constraints (such as observations of pulsar kicks, 
which constrain the supernova kick distribution) can be 
added by other means, as prior distributions on the space 
of model parameters. We fully expect to have much 
stronger constraints on our understanding of population 
synthesis in the near future. 

Stronger constraints, however, will require a consider- 
ably more systematic approach than the straightforward 
presentation we have used here. As the expected uncer- 
tainties decrease, greater care must be taken to include 
every uncertainty, no matter how minor, many of which 
for clarity we have neglected. For example, in future cal- 
culations we expect to include uncertainties in L m i n and 
p, our fit and rate estimates, and even the star forma- 
tion rate we use to convert simulation results into rate 
estimates. Additionally, we will self-consistently choose 
the constraint confidence intervals in order to construct a 
meaningful posterior confidence interval on each merger 
rate. 

Note on BH-pulsar rates: Recently. iPfahl et alJ J2005) 
have estimated that BH-millisecond pulsar systems 
should be exceedingly rare: they find an upper bound 
of 10~ 7 /year on the formation rate of binary systems 
containing a BH and a recycled pulsar. In our own 
computations, we have never seen such a system form, 
even though we have seen i=s 6 x 10 4 merging NS-NS 
binaries form. This result suggests the branching ra- 
tio for BH-PSR:NS-NS formation is < 10~ 4 . If we 
use a conservative value for the merging NS-NS forma- 
tion rate, 10~ 4 /yr/galaxy, then we expect the forma- 
tion rate of BH-PSR binaries to be significantly less than 
10~ 8 /yr/galaxy, entirely consistent with their constraint. 
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APPENDIX 

A. CALCULATING DCO EVENT RATES WITH POPULATION SYNTHESIS 

This paper relies upon formation rates extracted from a large sequence of archived population synthesis calculations. 
This appendix describes how these archives were generated and used to produce formation rates. It also explains the 
expected uncertainties in each merger rate estimate. 

How rates were estimated 

Principles of archive generation: For a given combination of population synthesis parameters, we generate a large 
collection N of binaries. To add a binary to the archive, we first generate progenitor binary parameters (mi, 7712, a, e) 
[i.e., the two progenitor masses (m^j) and the initial semimajor axis (a) and eccentricity (e)] which are (i) drawn from 
the distribution functions presented in OKB and which (ii) satisfy any conditions we impose to reject binaries irrelevant 
to the study at hand - e.g., mi,2 > 4, or more elaborate conditions described in OKB. (The latter conditions offer 
a significant speed improvement, but rely upon the experience gained in previous runs to insure that the conditions 
imposed do not reject physically relevant systems.) Given satisfactory initial conditions, binaries are assigned a 
randomly-chosen formation time, then evolved from whenever they form until the present day. Binaries are successively 
added to the archive until some termination condition is reached - typically, that the number n of a given class of 
binaries, such as merging NS-NS binaries, has crossed a threshold (i.e., n = 10) . 

Archive classes: Our population synthesis runs are summarized in Table 1X21 Each run of the population synthesis 
code falls into a certain class, depending on what choices were made for (i) the target sys tems on which the termination 
threshold was set (i.e., stop when we get n merging BH-BH binaries; see column 2 of Table lX2|) . (ii) the specific threshold 
n chosen (i.e., which insures that the formation rate of the target system type is determined to an accuracy roughly 
~ 1/y/n; see column 3 of Table IX2|) . and (iii) the combination of conditions applied to filter progenitor binary systems 
(see the last column of Table IX2|) . In Table IX2l the filters B, and S correspond to using the partitions presented in 
listed in Table 2 of OKB for BH-BH binaries (B) and NS-NS binaries (S); the 'W filter uses only the first NS-NS 
partition listed in Table 2 of OKB, which filters out WD progenitors. Note the first column of Table IA2I merely 
provides a label for the archive class. 

Applying archives: ^From the ratio of the number of binaries of a given type seen to the number of binaries in a run, 
modulo a normalization factor presented in OKB, we can calculate the formation rates for any binary type of interest. 
However, to avoid extreme biases associated with a poor choice of filter or stopping condition, we use only certain 
archives to estimate merger rates, as described in Table IA3I [In this table, all rates are total merger rates, with the 
exception of the last two rows, which correspond to the visible merging (v) and visible wide (vw) NS-NS binaries.] 

Understanding errors in rate estimates 

Tables IX2l and 1X51 provide the information needed to understand errors in our formation rate estimates. 

Example: The BH-BH formation rate estimate is produced from a single archive (b). Archive b is a collection of runs 
which stop when 10 merging BH-BH binaries are found; while the filters in this archive can prevent the formation of 
binaries involving NS, they do not significantly limit BH-BH binary formation. Therefore, archive b is ideally suited 
to estimate the BH-BH merger rate to an accuracy of order l/^/lO ss 30%. 



TABLE A2 

Classes of runs 
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TABLE A3 














Classes used for specific rates 
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b 




b' 
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BH-BH 
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BH-NS 
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NS-NS 
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X 
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NS-NS(v) 
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NS-NS (vw) 
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Example: The NS-NS formation rate estimate is produced from an amalgam of archives (a', a", 6', and c). None 
of these archives applies filters which prevent NS-NS formation, though one (a') applies filters which prevent the 
formation of nearly anything else. However, these archives do involve different termination criteria: the first two 
terminate when a large number of NS-NS binaries have formed, whereas the last two terminate when only 10 BH-BH 
or BH-NS archives have formed. If only the first two were used, we could guarantee the NS-NS merger rate to be 
known to within 10% accuracy. However, to augment our statistics, we additionally included binaries from b' and c; 
while these archives should usually have many NS-NS binaries (cf. Fig.|3|), we cannot guarantee any minimum number 
a priori. To simplify error estimates, we selected only those elements of b' and c with more than 30 merging NS-NS 
binaries. Thus, we expect our NS-NS rate to be known to within 20% accuracy. 

Example: The estimate for the visible wide NS-NS formation rate is the least accurate and most challenging calcu- 
lation we performed. Since visible wide NS-NS systems were rare and since we did not (unlike the BH-BH case) have 
a simulation dedicated to discovering them, we had to scavenge through all the archives which could have produced 
them (i.e., not b) in sufficient numbers (i.e., not a) to permit a moderately accurate rate estimate. In practice, we 
selected those runs which formed more than 8 wide visible NS-NS binaries, which should in principle give us an ac- 
curacy of order (few) x 35%. [In practice, we found a roughly 75% accuracy, comparable to the 65% accuracy of our 
next-most- accurate estimate (the BH-NS rate).] 

B. SAMPLE FITS TO MERGER RATES 

This paper and in OKB rely upon fits to 7-dimensional functions obtained from the StarTrack population synthesis 
code. In this section, we provide an example of an explicit formula for a quadratic- order polynomial fit to the BH-BH 
merger rate. We express this fit in terms of the following dimensionless parameters x^ € [0,1]: x\ = r/3 where r € [0, 3] 
is a negative power-law index describing the mass ratio distribution assumed for binaries; X2 = w characterizes the 
strength of stellar winds; the kick velocity distribution consists of two Maxwellian distributions with 1-D dispersions 
of o\ and 02, which are varied within [0,200kms _1 ] and (200, 1000kms _1 ], respectively, using two parameters: 
£3 = ci/(200km/s) and X4 = —1/4 + er 2 /800km/s); a third parameter X5 = s is used as the relative weight between 
low and high kick magnitudes; xq = aX is the effective common-envelope efficiency; and x-j — f a is the fraction of 

mass accreted by the accretor in phases of non-conserv ative mass transf er. [These parameters are discussed more 

thoroughly in OKB and in the original StarTrack paper IBelczvnski et alJ ((2002).] In terms of these parameters, we 
find the following quadratic fit to the BH-BH merger rate: 

log 10 [^BH-BHyr] = -5.84517 + 1.30448 x x ~ 0.406066 xi 2 

-0.310686 x 2 - 0.407175 xi x 2 + 0.0142072 x 2 2 
+0.717803 x 3 - 0.367487 x x x 3 - 0.48743 x 2 x 3 
+0.29931 x 3 2 - 0.7701 74 x 4 + 0.242792 x x x 4 
-0.0811259 x 2 x A - 0.0954582 x 3 x 4 + 0.113668 x 4 2 
+0.460929 x 5 - 0.114934xi x 5 + 0.490873 x 2 x 5 
-0.691954 x 3 x 5 + 0.588787 x 4 x 5 - 0.810968 x 5 2 
-3.27367 x 6 - 0.214978 x x x 6 + 0.674502 x 2 x 6 
-0.779896 x 3 x 6 + 0.0891919 x 4 x 6 + 1.37719 x 5 x 6 
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Fig. C7. — Cumulative probability distributions P^{X) denned so P^{X) is the fraction of all models consistent with both constraints 
with x k < X. The left panel shows the distributions for the 3 kick-related parameters 23 , X4 , £5 ; in this panel, the bottom (dashed) curve 
denotes P3, the middle (solid) curve denotes P4, and the top dotted curve denotes P5. The right panel shows the distributions for Pi (solid 
top line), P2 (dashed), Pq (dotted), and P7 (solid, bottom curve). 



+2.30296 x 6 2 + 1.68227 x 7 - 0.289592 x x x 7 
-0.19047x 2 x 7 - 1.18196x 3 x 7 - 0.0281177 x 4 x 7 

+0.5 1 70 42 x 5 x 7 + 0.6 7 5 96 x 6 x 7 - 0.45 4 39 x 7 2 (Bl) 

Relation of this fit to those used in paper: The fit presented above is substantially less accurate than those actually 
used in our paper or in OKB: it is accurate only to within a factor 1.8 ±:L (i.e., when we evaluate this fit at all of our trial 

1/2 

points, we find the standard deviation between our results and the fit to be ((loglZBH-BHjit — ^oE^-bh-bh) 2 ) = 
0.26). The fits actually used in this paper are typically cubic (120 parameters) and quartic (330 parameters) order. 
Given the large number of these parameters we chose to just provide the quadratic-order fit as an example above. 
However, the authors are happy to provide the much longer expressions for the higher-order fits used upon request by 
any reader. 

The rate functions are demonstrably not separable: we cannot fit the rate functions well with a function of form 
Xl(xi)X2(x2) . . .X 7 {x 7 ). Even this toy fit contains strong off-diagonal terms. 

C. CHARACTERIZING THE CONSISTENT REGION OF POPULATION SYNTHESIS MODELS 

Using our monte-carlo method to select models compatible with our constraints, we have found a relatively small 
seven-dimensional volume consistent with observations that corresponds to about 2% of all the runs we performed. 
Unfortunately - with some exceptions - this volumetric constraint does not translate to easily-understood and strong 
constraints on the individual parameters Xk (using the notation of Appendix |B|) . On the one hand, because the 
dimension is high, weak constraints on each parameter can correspond to very strong volumetric constraints. On the 
other hand, as demonstrated in Figure fHt 7 . because the consistent region is extended through our high-dimensional 
model space in a inhomogeneous anisotropic fashion, wide ranges of values of each parameter are still allowed, even 
after applying the constraints. 

To provide the reader with a global view of the parameter values associated with the models that turn out to be 
consistent with our constraints, in Figure IClf we show the cumulative distributions of the consistent model parameter 
values for each of the seven parameters. It is evident that the full ranges of values ([0,1]) are covered by the model 
parameters x^ for the set of models consistent with the constraints. However, certain qualitative conclusions can be 
drawn: about 80% of the consistent models have kick relative weights (parameter X5) below 0.3; about 50% of the 
consistent models have mass-ratio power-law indices smaller (in absolute value) than 0.6 (xi < 0.2; fractions of mass 
lost from the binary during non-conservative mass transfer phases in the range 20%-60% are not favored. 

Given the above, any simple attempt to describe the consistent region will necessarily be a crude approximation. 
Nonetheless, for completeness we attempt to characterize the extended consistent region through its mean values. The 
mean model consistent with our constraints is given by: 

x = (xi,x 2 ,...,x 7 ) = (0.33, 0.46, 0.61, 0.53, 0.18, 0.43, 0.57) (CI) 

These mean values correspond to a model with: fairly flat mass ratio distribution (power-law of ~ — 1); moderate 
stellar winds (strengths reduced by factors of ~ 2); moderate kicks drawn from Mawellians with ci ~ 120 km s -1 , 
(T2 ~ 625 km s -1 , and with relative weights of ~ 20% (favoring u% over <7i); moderate values for an effective common- 
envelope efficiency (aX ~ 0.4 including the central concentration parameter A); and moderately non-conservative mass 
transfer phases (~ 40% of the mass is lost from the binary). It is very important though to keep in mind the broad 
ranges of these parameters shown in Figure IClf. 



